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opendf - an implementation of the dual fermion method for strongly correlated systems 

Andrey E. Antipov, 1 ’* James P.F. LeBlanc, 1 and Emanuel Gull 1 
1 Department of Physics, University of Michigan, Ann Arbor, Michigan 48109, USA 

The dual fermion method is a multiscale approach for solving lattice problems of interacting 
strongly correlated systems. In this paper, we present the opendf code, an open-source implementa¬ 
tion of the dual fermion method applicable to fermionic single-orbital lattice models in dimensions 
D = 1,2,3 and 4. The method is built on a dynamical mean field starting point, which neglects 
all local correlations, and perturbatively adds spatial correlations. Our code is distributed as an 
open-source package under the GNU public license version 2. 


I. INTRODUCTION 

Understanding the physics of complex correlated electron systems beyond simple approximations or exactly solvable 
limits is a long-standing goal of condensed matter physics. Towards this goal, the dynamical mean field theory 
(DMFT) 1 6 is a workhorse which provides numerically simulated results for the physics of such systems. It establishes 
that, if correlations and interactions are assumed to be local, the (intractable) extended system can be mapped 
self-consistently onto an Anderson impurity model, which can then be solved numerically. 

The dynamical mean field approximation of locality is often precise enough that general material trends can be 
reproduced. Nevertheless, cases where non-local correlations lead to behavior not captured by DMFT are known ~ 10 , 
and therefore methods that improve on this approximation are needed. The dual fermion method 1 , which perturba¬ 
tively adds corrections to a DMFT starting point reintroduces momentum dependent correlations. If all corrections 
are included, the method recovers the full momentum dependence of the original problem and becomes numerically 
exact. 

In this paper, we present opendf, an implementation of the ‘ladder series’ variant of the Dual Fermion method . 
This variant is approximate, as neither vertices with more than four legs nor series of vertices beyond a single ladder 
are considered. Nevertheless, it has been shown to consistently improve on DMFT results “ and capture critical 
properties of phase transitions 1 ’ . 

Dual fermion calculations rely on a dynamical mean field input, which can be provided by one of the publicly 
available open source software packages that implement the approximation, including ALPS (the Algorithms and 
Libraries for Physics Simulations) 1 ”, TRIQS (the Toolbox for Research on Interacting Quantum Systems) 2 ' 1 , and 
iQIST . This initial step requires the self-consistent solution of an interacting quantum many-body system and the 
calculation of vertex functions* ’ and is computationally much more expensive than the summation of the dual 
fermion diagrams. 

The rest of this paper is organized as follows: Section II introduces the methodology. Section III describes distri¬ 
bution aspects, section IV performance aspects, section V shows some examples, and section VI will conclude. 


II. METHODOLOGY 
A. Prerequisites 

We consider a general fermionic single-orbital lattice model with a Hamiltonian 
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written in mixed momentum, k , and real space, i, notation in terms of creation and annihilation operators (cl and 
Cka respectively). The index a labels the spin projection, e& is the lattice dispersion relation and k is the vector in 
the reciprocal space. H lnt is the local interaction for each site, i. on the lattice. No assumption is made within DF as 
to the structure of H lnt . 

As a first step, which must be performed outside of this code, an approximate solution of the model is obtained from 
a dynamical mean field calculation, for example provided by the ALPS code with an appropriate impurity solver . 
It provides an estimate for the local Green’s function of the lattice problem as a solution of the Anderson impurity 
model, embedded into a self-consistently determined hybridization. The imaginary time action of this “impurity 
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problem” reads 
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where = f/ dr_ff lnt [c] (r), Cj(r)] is the interaction part of the action and A W(T is a self-consistently determined 
hybridization function. The DMFT impurity solver computes the one particle Green’s function g u<7 = -(c^c^) of 
the of the Anderson impurity model and the two particle vertex functions (i.e. the connected parts of two-particle 
Green’s functions) 
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The following quantities are then provided as an input to the DF simulation: 

• g u - the full Green’s function of the DMFT impurity problem (same values for both spin components) 

• A u - hybridization function of the DMFT impurity problem 

• g - chemical potential of the problem 

• Two independent components of the impurity vertex function, , from Eqn (3): and 


TUt — U 

The present version of the code considers only spin-symmetric solutions of fermionic spin s = 1/2 problems, and 
does not describe symmetry-broken phases. We will omit the spin index a in single particle quantities in what follows. 


B. Ladder dual fermion self-consistency loop 


A precise derivation of the DF equations can be found in . Here we outline the equations solved within the 
opendf code. The evaluation of the DF equations starts with the construction of the bare dual fermion propagator 

= [9u 1 + - £fc] ^ gu, (4) 

which represents a /c-dependent correction to the impurity Green’s function. This Green’s function is used to construct 
two-particle bubbles: 
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Here the integral over the Brilloin zone is replaced with a discrete summation with Ni~ points in each direction. The 
impurity vertex functions are combined into density and magnetic channels (labeled d/m respectively) as: 
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The vertices for the respected channels from Eq. 6 and the bubbles from Eq. 5 are substituted into ladder equations: 
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r[j iU y is called the fully dressed vertex function. 

Evaluation of Eq. 7 is performed independently for each pair of bosonic frequencies fl and transfer momenta q. 
ln,u,u' and rn,^, w '(g) are represented as matrices in the space of fermionic Matsubara frequencies ui, a/, and \n,u"(q) 
is a diagonal matrix. In this matrix notation, Eq. 7 reads 

(i - 7x)f = 7. (8) 


This equation is physically correct only when the maximum eigenvalue of 7 % is smaller than one, i.e. all eigenvalues 
of the matrix D = 1 — 7 X are positive. Eq. 8 is then solved and T is obtained. When the determinant of D is negative 
and a negative eigenvalue exists, the DF solution is outside of the convergence radius of the ladder approximation. 
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Nevertheless, given that the resulting solution is unique, one can extend this convergence radius by doing a low-order 
iterative evaluation of T and checking if the inversion of Eq. 7 can be obtained on the next DF iteration. 

Once the fully dressed vertex function Tn iWiU / is obtained, it is used in the Schwinger-Dyson equation to obtain the 
dual self-energy E Wi fc. The equation reads: 
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where T^ 2 ) = 7X7 indicates the second order (first iteration) correction from Eq. 7 to avoid diagrammatic double 
counting. 

The resulting dual self-energy is used to obtain the dual Green’s function from the Dyson equation: 
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The procedure is repeated until convergence of G is achieved. 


C. Resulting observables 

The fully converged dual Green’s function G, self-energy E, vertices T d / m determine the lattice correlators. Specif¬ 
ically, 

• the lattice self-energy: 
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where E FMFT = ito + /z - A w - g u 1 . 

• The lattice Green’s function 

= [A w - £ fe ] _1 + [A w - \~ x [A u - e fc ] _1 . 

Eqs. 11 and 12 are related by a Dyson equation for G and E. 

• The charge and spin susceptibilities 
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III. DISTRIBUTION 
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The dual fermion code is distributed as a C++ library with compiled executables hub_df _cubicDd, where D labels 
the number of dimensions (D = 1,2, 3,4). We use the opensource gftools library for algebraic operations with 
single- and multi-particle Green’s functions and its interface to the ALPSCore libraries' for loading/saving hdf5 
objects. The code and the documentation are available as Ref.' . 
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IV. EXAMPLE I AND PERFORMANCE ANALYSIS 


As a first example and illustration of the performance of the code, we provide an example input generator for the 
particle-hole symmetric Hubbard model at U 3> t (the “atomic limit”) . In this case the input quantities are given 
analytically by: 
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where A u = 1 + /7 2 /(4w 2 ) and uq = cquq = w + fl,W 3 = u/ + f 2 ,W 4 = u/ is used to simplify the notation. The 
corresponding program is provided with the code. 

The numerical solution of dual fermion equations requires introducing several control parameters. In particular, the 
vertex function 7 n ,u,u' is sampled on a grid with a cutoff An in bosonic and N u fermionic frequencies and the Brilloin 
zone is sampled on a finite grid of size Nk, giving a total volume of the system of Njp. We analyze the convergence 
of the code upon tuning TVn, N u and Nk and the computational effort below. Eqs. 16, 17, 18, 19 are used to provide 
the input to the code and the system is evaluated in 2 dimensions, at U = 20, /j, = P = 1 . We choose the value of 
g = G in /p, 0,0 t° control the convergence. We then plot the normalized difference 
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as a function of control parameter N Xl with x = {fi, w, fc}, and extrapolate N x —> 00 to evaluate the error. For the 
most expensive point shown here, the run-time of the simulation was ~ 2 min on a laptop. 


2d, U=20, (3 = 1,48 fermionic freqs, 16 2 kpoints 



FIG. 1. (a) Execution time of the dual fermion calculation for the Hubbard model in 2 dimensions with “atomic limit” input 
at U = 20, P = 1 as a function of the number of bosonic frequencies An at = 48, Nk = 16; (b) Systematic error S g of the 
dual fermion Green’s function Giui,k at ioj = wr/ p, k = (0, 0) as a function of bosonic frequencies An, plotted on a logarithmic 
scale. 


Fig. 1 shows the performance of the opendf code upon the change of the total number of bosonic frequencies Aq 
in the vertex 70 for a fixed number of fermionic frequencies A w = 48 for a 16 x 16 k-space grid. The computational 
effort, indicated by the time to convergence in Fig. 1(a), grows linearly in Nq. The error 5 g , as defined in Eqn. 20 
and shown in frame (b), is of the order of a percent and decreases with a power law. 
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2d, U=20, (3 = 1,3 bosonic freqs, 8 2 kpoints 




FIG. 2. (a) Execution time of the dual fermion calculation for the Hubbard model in 2 dimensions with “atomic limit” input 
at U = 20, p = 1 as a function of the number of fermionic frequencies at An = 3, Nk = 8; (b) Error 5 g of G; w ,fc at 
iui = in/P,k = (0, 0) as a function of 1 /N u , for the same parameters. 


We analyze the performance of the code with respect to the change of the total number of fermionic frequencies N u 
in Fig. 2. In this benchmark we fix the number of bosonic frequencies, = 3, and perform the calculation on a 8 x 8 
k-space grid. The computational expense seen in Fig. 2(a) grows almost quadratically, while the relative error shown 
in Fig. 2(b) is an order of magnitude smaller, as compared to the variation in Nq shown in Fig. 1(b) and reduces as a 
power-law with an increase of N u . 


2d, 11=20, (3 = 1,3 bosonic freqs, 48 fermionic freqs 
16 2 32 2 48 2 64 2 




FIG. 3. (a) Execution time of the dual fermion calculation for the Hubbard model in 2 dimensions with “atomic limit” input 
at U = 20, P = 1 as a function volume N/P at Nn = 3, N u = 48. (b) Systematic error 5 g of Giu,,k at iui = in/P, k = (0, 0) as a 
function of the volume at the same parameters. 


The performance of the code with respect to the change of number of k-space samples within the Brilloin zone Njp , 
is plotted on Fig. 3. The computational effort (frame (a)) scales linearly with the volume Njp of the system and 
shows fast convergence of the relative error 5 g (frame (b)). 


V. EXAMPLE II - HUBBARD MODEL, 2 DIMENSIONS 

We provide the practical illustration of the method for the Hubbard model in D = 2 dimensions. We show the 
/c-dependence of the real part of the lattice self-energy E(fc, *w n ) at ito n = iuio = in/p in Fig. 4 for the case of 
particle-hole symmetry at U/t = 8 and compare it with available data from the Dynamical Cluster Approximation . 
The impurity model, solved using the ALPS DMFT package with a CT-AUX solver ! , was used as an input. The 
DMFT self-energy is momentum-independent, ReX^ IFT = 0, and is plotted with a dashed line. Taking into account 
the spatially dependent corrections by the dual fermions leads to a correct momentum-dependence of the self-energy, 
matching in this case the DCA result. A detailed comparison between multiple methods will be discussed elsewhere . 
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FIG. 4. Momentum dependence of the real part of the lower Matsubara frequency of the lattice self-energy of the particle- 
hole symmetric Hubbard model in 2 dimensions, as obtained by the opendf calculation at U/t = 8, j3 = 0.5 along the 
( k x ,k y ) = (0,0) —» (7r, 0) —> (n,7 r) —» (0,0) path (red points). Shown also is the DMFT value (as a dashed line) and the 
comparison data from the 72-sitc Dynamical Cluster Approximation calculation (solid blue line). 


X ch (Q = 0) —— X sp (^ = 0) 



FIG. 5. X ch (left panel) and y sp (right), the static spin and charge susceptibilities, at fi = 0 as a function of momenta q x and 
q v as obtained by the opendf calculation for /? = 0.5, U/t = 8. 


We illustrate the susceptibility in Fig. 5. Plotted are the static spin- and charge- susceptibilities at U/t = 8 for the 
particle-hole symmetric case. The spin susceptibility, peaked at (w, ir) due to antiferromagnetic fluctuations is much 
larger than the charge one. 


VI. CONCLUSION 


In this paper we have introduced an open source implementation of the dual fermion method, the opendf project. 
It solves the dual fermion self consistency equations and computes non-local corrections to the local solutions provided 
by DMFT. opendf can be used to augment DMFT computations with two-particle quantities and add momentum 
dependence to DMFT observables. 

Future development of the code is anticipated. Further releases will include extensions to additional diagrams, 
broken-symmetry phases and multi-orbital systems. 
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